STEM-EDXS spectrum image analysis with HyperSpy¶
Supplementary information to the paper: "Structural and chemical properties of superconducting Co-doped BaFe$_2$As$_2$ thin films grown on CaF$_2$"
This notebook was run with the HyperSpy (1.5.2) python package. For information on HyperSpy you can visit the website: https://hyperspy.org/
Information about the usage of Jupyter notebooks and how to run this notebook interactively: https://jupyter.org/. Put the ".ipynb" file and the data "STEM-EDXS-data_Fig4b.bcf" in the same folder.
This notebook shows the data treatment for Figure 4b (sample/slow). The other shown STEM-EDXS datasets were treated in a similar way. The data was collected on an FEI Tecnai Osiris equipped with ChemiSTEM technology with Bruker's Esprit software. The EDXS spectrum image is saved as a Bruker composite file (bcf) which can be imported by HyperSpy. Principal component analysis (PCA) is used to reduce noise. Qualitative elemental maps, i.e. background-corrected peak intensities at each pixel, are extracted from the denoised datacube.
Preparation¶
We load the necessary packages. Here the "matplotlib nbagg" GUI backend is used, but you can also use the "qt" backend, which creates pop-up windows.
%matplotlib nbagg
#%matplotlib qt
import hyperspy.api as hs
import numpy as np
import matplotlib
import os
Create additional subdirectories to save results:
def create_folder(savepath):
try:
os.mkdir(savepath)
except OSError:
print ("Creation of the directory %s failed" % savepath)
else:
print ("Successfully created the directory %s" % savepath)
#Main directory for results
HSpath = 'HyperSpy Analysis'
#Elemental maps subfolders
PCApath1 = HSpath + '\\Elemental maps'
PCApath2 = PCApath1 + '\\png'
create_folder(HSpath); create_folder(PCApath1); create_folder(PCApath2);
Load the datafile:
file = "STEM-EDXS-data_Fig4b.bcf" #(bcf) filename
s = hs.load(file, select_type='spectrum_image')
First we inspect the original data. Note that the spectrum image was acquired rotated by 90° compared to the figure in the paper: On the left is the CaF$_2$ region. In the middle is the Ba122 film. The bright strip is the Pt protective layer.
s.plot() #Inspect data
Save summed-up spectrum and metadata of the STEM-EDXS spectrum image:
s.sum().save(HSpath+'\\'+'SumSpectrum.msa') #mean spectrum
s.original_metadata.export(HSpath+'\\'+'Metadata.txt') #metadata
Data pretreament before PCA¶
PCA is sensitive to outliers/periodic noise which can be falsely interpreted as a "real" signal. Data treatment (rebinning of pixels, outlier removal, ...) prior to PCA is therefore recommended. In this case, we crop away the internal detector peak at 0 keV and unnecessary energies above 16.5 keV (no X-ray line peaks).
s.sum().plot() #Plot sum spectrum to see relevant energy range
startE = '0.2 keV'
endE = '16.5 keV'
s = s.isig[startE:endE]
Cropped energy range:
s.sum().plot() #Plot sum spectrum
Multivariate statistical analysis¶
In the next steps we perform PCA to enhance the signal-to-noise ratio of our dataset.
Perform PCA decomposition¶
For the PCA calculation we convert the integer values (X-ray counts) to float values:
s.change_dtype('float') #Convert integer (X-ray counts) to float
Now we run PCA. We assume Poissonian noise in the EDXS dataset. For the decomposition, the singular value decomposition (SVD) algorithm is used.
s.decomposition(normalize_poissonian_noise=True)
After the decomposition has finished, we visualize the results by running the next cell. 3 windows will open:
- The "scree plot" is used to determine the number of relevant components. See also, e.g. https://en.wikipedia.org/wiki/Scree_plot. A vertical line will be close to the elbow position, which is an assumption made by HyperSpy. However, this may not always be correct. Therefore we also use the score maps/loadings to find the cut-off point for non-noise PCA components.
- The score maps and (abstract) component spectra are in the other two windows: When you click on them you can use the arrow keys to scroll through the individual score maps/spectra.
Here, the maps show features for index 0 to 4 and only noise starting from 5. So in this case the guess by the vertical line is correct and we want to keep the first 5 components.
ax = s.plot_explained_variance_ratio(vline=True) #Close windows after inspection
ax2 = s.plot_decomposition_results() #Close windows after inspection
Specify number of meaningful components in the dataset¶
By inspection of the scree plot/score maps we can decide on the number of components (here 5):
N = 5 #Specify output dimensions/number of meaningful components
Reconstruct datacube with $N$ meaningful components¶
PCA_filtered = s.get_decomposition_model(N)
We can visualize original data (red) and noise-reduced data (blue) at each pixel. The PCA-filtered spectrum is greatly enhanced compared to the raw X-ray counts (on a single pixel) in the original dataset.
(s + PCA_filtered * 1j).plot()
Extract elemental maps¶
In this part, we continue by extracting qualitative maps from the PCA-filtered datacube. We start by saving the calibrated dimensions of the EDXS mapping to a text file:
#Read spatial calibration from metadata
PixelsizeInNm = s.original_metadata.Microscope.DX*1e3 #from bruker bcf, metadata is in µm -> *1e3 for nm
Pixels = [np.shape(PCA_filtered)[0], np.shape(PCA_filtered)[1]] #pixels in x,y
MapsizeInNM = [PixelsizeInNm*np.shape(PCA_filtered)[0], PixelsizeInNm*np.shape(PCA_filtered)[1]] #total mapsize in nm
#save to a txt file
f = open(PCApath1+"\\Scale.txt", "w")
f.write("Scale for qualitative EDS mappings:\n")
f.write("Dimensions:\t\t %.0f" % np.shape(PCA_filtered)[0] + " x %.0f" % np.shape(PCA_filtered)[1] + "\n")
f.write("Pixelsize in nm:\t %.3f" % PixelsizeInNm + '\n')
f.write("Total size:\t\t %.2f" % MapsizeInNM[0] + " x %.2f" % MapsizeInNM[1] + " nm\n")
f.close()
Case 1: No X-ray line peak overlap¶
We define background and integration windows for all X-ray lines of interest and then extract the intensities. Background-window (and integration-window) positions are defined for each X-ray line before extracting the background-subtracted intensities. The variable "iww" is the width of the integration window for signal extraction. "xraylines" are the desired X-ray lines to extract.
iww = 1.2 # width of integration window times FWHM of X-ray peak of interest
xraylines = ['As_Ka', 'Ba_La', 'Ca_Ka', 'Fe_Ka', 'O_Ka'] #Sort alphabetically
Now we estimate the background/integration windows to get their rough placement on the energy axis:
#Dont change anything here
dummy = xraylines[0].split('_') #Get element name without X-ray line, e.g. "As" from "As_Ka" in e
PCA_filtered.set_elements([dummy[0]]) #Dummy line to erase all elements already present
PCA_filtered.set_lines(xraylines) #Set desired X-ray lines
bw = PCA_filtered.estimate_background_windows(line_width=[1.5, 1.5], windows_width=0.6)
iw = PCA_filtered.estimate_integration_windows(windows_width=iww)
We look at the summed-up spectrum and plot the current background/integration windows. Leave the plot open, we will use it to fine-tune the window positions.
PCA_filtered.sum().plot(True, background_windows=bw, integration_windows=iw)
We use the "Zoom" and "Pan" tools in the figure to inspect the X-ray lines. For each X-ray line of interest, you should see (a) the two background windows (filled lines), (b) the integration windows (dashed lines), and (c) a black line indicating the currently assumed background. The background is modeled as a straight line which runs through the mean values of the two outer background intervals. The "bw" variable with the background window position is an array with a row for each X-ray line with 4 entries: left/right coordinate of left window and left/right coordinate of right window (in keV):
bw
In this example we update the O-K$\alpha$, As-K$\alpha$ and the Ba-L$\alpha$ background window positions. For example, here are the old positions for O-K$\alpha$:
PCA_filtered.sum().plot(True, background_windows=bw, integration_windows=iw) #Pan to O-Ka line
Now we define our new background window positions and overwrite the values:
#Define new background positions here
bw_AsKa = [10.17 , 10.3 , 12 , 12.3]
bw_BaLa = [4.17, 4.28, 6.01, 6.18]
bw_OKa = [0.4, 0.42, 0.86, 0.89]
#Overwrite old positions with new ones
bw[0,:] = bw_AsKa
bw[1,:] = bw_BaLa
bw[4,:] = bw_OKa
Plot again to check updated positions:
PCA_filtered.sum().plot(True, background_windows=bw, integration_windows=iw) #Pan to O-Ka line
The O-K$\alpha$ line partially overlaps with other X-ray lines in the low-energy region of the EDXS spectrum but the extraction by this window method gives a reasonable O elemental map. Later in this notebook, we will compare it to the map obtained by peak-fitting the O-K$\alpha$ line.
Extract intensity (dashed line region) at each pixel:
intensity = PCA_filtered.get_lines_intensity(background_windows=bw, integration_windows=iw)
intensity
For example, we can inspect the qualitative map for As-K$\alpha$ (index 0):
intensity[0].plot()
Save maps to files:
for i in range(np.shape(bw)[0]):
np.savetxt(PCApath1+"\\"+str(PCA_filtered.metadata.Sample.xray_lines[i])+".txt", intensity[i].data[:,:])
matplotlib.image.imsave(PCApath2+"\\"+str(PCA_filtered.metadata.Sample.xray_lines[i])+".png", intensity[i].data[:,:])
Case 2: Fit Gaussian functions to seperate overlapping X-ray lines¶
In case of peak overlap we can use Gaussian functions to model the X-ray line families by multiple linear least-squares fitting. This is used here for the separation of F K$\alpha$ (0.677 keV) and Fe L$\alpha$ (0.705 keV) X-ray lines. The following procedure is used:
- Crop energy dimension to overlap region to speed up computation.
- Initialize a HyperSpy model: Relevant X-ray lines in the energy region will be added automatically as Gaussian functions at their respective X-ray line energy.
- Add a background function to the model: Here a linear funtion is used (polynomial of 1st order).
- Fit the model at specific locations to check the fit and to initialize some sensible starting parameters.
- Run "multifit" to run the fitting procedure at each pixel.
- Extract the amplitude parameter of the relevant Gaussian functions to generate elemental maps.
Crop signal to critical overlap region. Here we also include the O K$\alpha$ line.
PCA_filtered = PCA_filtered.isig[0.4:0.9]
Plot the cropped energy range:
PCA_filtered.sum().plot()
Set up X-ray lines and fit function: "beam_energy" should be set to the upper limit of the (cropped) energy range (in keV), so that Hyperspy adds the correct X-ray line models. We then specify elements which are expected to be found in the sample. Spurious signal from the microscope enviroment, e.g. Cu or Zr, or other samples areas (e.g. Pt from deposited Pt, Ga from FIB preparation) are also added.
PCA_filtered.set_elements([]) #Delete X-ray lines/elements, which may be stored in the bcf file
PCA_filtered.set_microscope_parameters(beam_energy=0.9)
PCA_filtered.add_elements(['Al', 'As', 'Ba', 'C', 'Ca', 'Co', 'Cu', 'F', 'Fe', 'Ga', 'Si', 'Ni', 'O', 'Pt', 'Zn', 'Zr'])
Initialize the model: HyperSpy will automatically add all X-ray lines from the specified elements which are below beam_energy.
m = PCA_filtered.create_model(auto_background=False)
Add a background function:
m.add_polynomial_background(order=1) #Linear function of order 1 background for TEM
In this energy range the background seems to have a negative slope. We help the fitting procedure by restricting the fit range for the slope variable to negative values.
m.background_components[0].a1.bmax = 0 #Set maximum value of slope to 0, i.e. only negative slopes are allowed during the fit
Print all functions/parameters of the model, check here if all lines of interest are present.
m.print_current_values(fancy=False)
In the EDXS fit some assumptions are used for fitting: The positions ("centre") and widths ("sigma") of the Gaussians are fixed, which is indicated by the "False" keyword in the column "Free". Furthermore, the relationships between peak heights in a X-ray line family are fixed (line weights). As a result, the amplitude parameter ("A") of the strongest X-ray line in a family is fitted and the amplitudes of other X-ray lines in the same family are scaled according to the X-ray line weights.
Note: HyperSpy 1.5.2 is missing the Ba M-line family, which was added by hand as Ba-Ma, Ba-Mb, ... lines in the "...\Anaconda3\Lib\site-packages\hyperspy\misc\elements.py" file. The Siegbahn nomenclature for the X-ray line names is not correct here. The values for the X-ray line energies and line weights were taken from the EPQ library by NIST: https://github.com/usnistgov/EPQ/tree/master/src/gov/nist/microanalysis/EPQLibrary.
We check if the fit looks reasonable at different sample positions:
#Plot the model
m.plot(True)
We are at position (0,0), i.e. top left. This is the CaF$_2$ region. We can fit our model at this position:
m.axes_manager.indices = (0, 0) #Move to (0,0)
m.fit(bounded=True)
In CaF$_2$ the total fit looks good (blue line on red dots). The two larger Gaussian functions are from F-K$\alpha$ (purple) and O-K$\alpha$ (black). We also check fit in the Ba122 region:
m.plot(True)
m.axes_manager.indices = (20, 20)
m.fit(bounded=True)
The fit in Ba122 is not perfect, which may be a result of a missing spurious X-ray line or some other artefact. However, for the separation of F-K$\alpha$ and Fe-K$\alpha$ this is seems not to be problematic. In the Ba122 region the Fe-K$\alpha$ peak (yellow) is clearly visible. Next, we assign the fit parameters of the last fit to all pixels to help the fitting procedure.
#Assign starting parameters of currently fitted position to all pixels of the spectrum image
m.assign_current_values_to_all()
Now we run a multifit, i.e. fit on every pixel ($70\cdot 70 = 4900$ positions). Depending on the CPU this can take a few minutes. A progress bar will be shown.
m.multifit(fitter="leastsq", bounded=True) # May take a long time! A progress bar will be displayed.
After the fit has completed, we set up and extract intensities (amplitudes of the Gaussians) from our model. Adjust the X-ray lines in "e" accordingly. In this example, the intensities of F-K$\alpha$, Fe-L$\alpha$, and O-K$\alpha$ will be extracted.
e = ['F_Ka', 'Fe_La', 'O_Ka'] #Sort alphabetically
dummy = e[0].split('_') #Get element name without X-ray line, e.g. "As" from "As_Ka" in e
PCA_filtered.set_elements([dummy[0]]) #Erase present element list
PCA_filtered.set_lines(e)
mI = m.get_lines_intensity(e) #Get the A parameter, i.e. amplitude of the Gaussians of interest
Extract maps and save results. The filename has an added '-fit' to indicate is was extracted by model fitting.
for i in range(0,np.shape(e)[0]):
np.savetxt(PCApath1+"\\"+str(PCA_filtered.metadata.Sample.xray_lines[i])+"-fit.txt", mI[i].data[:,:])
matplotlib.image.imsave(PCApath2+"\\"+str(PCA_filtered.metadata.Sample.xray_lines[i])+"-fit.png", mI[i].data[:,:])
We can see if the fit went well by comparison of Fe-K$\alpha$ map extracted by the window method and the Fe-L$\alpha$ line obtained by fitting. Both Fe maps look similar, so the separation of Fe and F worked out fine. The faint F signal in the Pt region might result from scattering of electrons from the heavy Pt into the CaF$_2$ region, thereby exciting F X-ray lines.
hs.plot.plot_images([intensity[3], mI[1], mI[0]], label=['Fe-Ka (window)' , 'Fe-La (fit)', 'F-Ka (fit)'],
cmap='viridis', axes_decor='off',
colorbar='multi', scalebar=[0],
scalebar_color='white', suptitle_fontsize=16)
The peak-fitting procedure also produced a less noisy O map. The right vertical strip is the oxidized surface of Ba122. The more faint left vertical line is the interface of CaF$_2$ and Ba122.
hs.plot.plot_images([intensity[4], mI[2]], label=['O-Ka (window)' , 'O-Ka (fit)'],
cmap='cividis', axes_decor='off',
colorbar='multi', scalebar=[0],
scalebar_color='white', suptitle_fontsize=16)
